library( "igraph" )
library( "vegan" )
library( "e1071" )

##
## Standard one-mode projection
##
## Similarity value is the size of the intersection of two nodes' neighbor sets.
## E.g. if node A is has neighbors A:{a, b, c}, and B:{b, c, d}, their
## similarity is 2 (b & c).
##
## The value is normalized to 0-1 by dividing over the maximum similarity in
## the network.
##
omp.standard <- function( adj ) {
	w <- adj %*% t( adj )
	## normalize
	w <- w / max( w )
	
	return( w )
}

##
## Jaccard Index
##
## Similarity is the size of the intersection of two nodes' neighbor lists over
## the size of the union of the lists.
## E.g. for nodes A:{a, b, c} and B:{b, c, d}, similarity is 2/4 = .5.
##
omp.jaccard <- function( adj ) {
	w <- vegdist( adj, method="jaccard", upper=TRUE, diag=TRUE )
	w <- as.matrix( w )
	## covert to similarity
	w <- 1 - w
	
	return( w )
}

##
## Pearson's correlation
##
## Similarity between two items is relative to the mean of all other pair-wise
## similarities (?).
##
omp.pearson <- function( adj, cutoff=0 ) {
	adj <- t( adj )
	w <- cor( adj, method="pearson" )
	
	## remove edges less than cutoff
	w[ w <= cutoff ] <- 0
	
	return( w )
}

## Same as Jaccard (?)
omp.binary <- function( adj, cutoff=0 ) {
	w <- dist( adj, method = "binary", diag=TRUE, upper=TRUE )
	w <- as.matrix( w )
	w <- 1 - w
	
	return( w )
}

##
## Hamming distance (as similarity)
##
## Similarity is the number of changes needed to make two adjacency vectors the
## same.  E.g. for possible neighbors {a, b, c, d, e} and nodes A:{a, b} and
## B:{b, c, d} with adjacency vectors A:{1, 1, 0, 0, 0} and B:{0, 1, 1, 1, 0},
## their Hamming distance would be 3 since there are 3 differences.
##
## The value is normalized to 0-1 by dividing over the maximum similarity in
## the network.
##
omp.hamming <- function( adj ) {
	w <- hamming.distance( adj )
	## normalize
	w <- w / max( w )
	## covert to similarity
	w <- 1 - w
	
	return( w )
}

##
## Walktrap community
##
comm.walktrap <- function( g ) {
	comm <- walktrap.community( g, weights=E( g )$weight, modularity=TRUE, merges=TRUE, labels=TRUE )
	memb <- community.to.membership( g, comm$merges, steps=which.max( comm$modularity ) )
	
	## plot dendrogram
	dev.new()
	dend <- as.dendrogram( comm, use.modularity=TRUE )
	plot( dend, main="Walktrap", nodePar=list(pch=c(NA, 20)) )
	
	return( memb )
}

##
## Fastgreedy community
##
comm.fastgreedy <- function( g ) {
	comm <- fastgreedy.community( g, weights=E( g )$weight, modularity=TRUE, merges=TRUE )
	memb <- community.to.membership( g, comm$merges, steps=which.max( comm$modularity ) )
	
	return( memb )
}

##
## Spinglass community
##
comm.spinglass <- function( g, spins=2 ) {
	comm <- spinglass.community( g, weights=E( g )$weight, spins=spins )
	
	return( list( membership=comm$membership ) )
}

##
## Betweenness community
##
comm.betweenness <- function( g ) {
	comm <- edge.betweenness.community( g )
	memb <- community.to.membership( g, comm$merges, steps=12 )
	
	return( memb )
}

##
## Leading eigenvector community
##
comm.eigenvector <- function( g ) {
	comm <- leading.eigenvector.community.step( g )
	memb <- community.le.to.membership( comm$merges, steps=1, comm$membership ) 
	
	return( memb )
}

## Read in the file.
adjlist <- read.table( "C:\\Documents and Settings\\saxman\\workspace\\BioStats\\network.txt", sep="\t", header=TRUE )
g <- graph.data.frame( adjlist, directed=FALSE )

## Collapse multiple edges into a single weighted edge.
## However, get.adjacency doesn't seem to assign weights in the matrix.
E( g )$weight <- count.multiple( g )
## Remove multiple edges.
g <- simplify( g )

## Classify each vertext as either a patient or a drug.
V( g )$type <- is.bipartite( g )$type
## Build two one-mode projections, a patient network and a drug network.
proj <- bipartite.projection( g )

## Assign names to the vertices in the projections.
V( proj[[1]] )$name <- V( g )$name[ !V( g )$type ]
V( proj[[2]] )$name <- V( g )$name[ V( g )$type ]

## Assign colors to vertices.
V(g)$color <- ifelse( V(g)$type, "blue", "red" )

## Draw the bi-partite, partient-drug network.
lay <- layout.fruchterman.reingold( g, weights=E( g )$weight )
dev.new()
plot( g, layout=lay, vertex.size=3.5, vertex.label=V( g )$name, asp=FALSE, edge.color="gray", edge.width=1 )

## Create an NxM, drug-by-patient, network, using the projections to determine the axes.
adj <- get.adjacency( g )
adj <- adj[ rownames( get.adjacency( proj[[2]] ) ), rownames( get.adjacency( proj[[1]] ) ) ]
## if we know the names, we could use those to constrain the axes.
#adj <- adj[ -grep( "^RA", rownames( adj ) ), grep( "^RA", colnames( adj ) ) ]

## Distance/similarity (as similarity).
## Useful functions for testing are standard, jaccard, and hamming.
## pearson is in development.
w <- omp.jaccard( adj )

g <- graph.adjacency( w, mode="undirected", weighted=TRUE, diag=FALSE )
## could also use a complete matrix/graph and remove unwanted vertices.
#g <- delete.vertices( g, grep( "^RA", V( g )$name ) )

## Assign degrees and remove detached vertices.
V(g)$degree <- degree( g )
g <- delete.vertices( g, V(g)[ V(g)$degree == 0 ] )

## Community finding algorithm.
## Useful functions for testing are walktrap, fastgreedy, and spinglass.
## betweenness and eigenvector are in development.
memb <- comm.betweenness( g )

## Plot one-mode network.
colbar <- rainbow( max( memb$membership ) + 1 )
col <- colbar[ memb$membership + 1 ]
col[ is.na( col ) ] <- "grey"

lay <- layout.fruchterman.reingold( g, weights=E( g )$weight )
dev.new()
plot( g, layout=lay, vertex.size=3.5, vertex.label=V( g )$name, vertex.color=col, asp=FALSE, edge.color="gray", edge.width=1 )
# edge.label=E( g )$weight

## Hierarchical clustering.
dist <- as.dist( 1 - w )
c <- hclust( dist, method="ward" ) #members=memb$membership # could prime hclust with output of comm algorithm.
dev.new()
plot( c )

## Create five clusters from the dendrogram.
cm <- cutree( c, k=5 )

## Convert hclust clusters to a membership vector.
memb <- list( membership=c() )
for ( i in 1:length( V(g) ) ) {
	memb$membership[i] <- cm[[i]] - 1
}

## Colorize the hclust clusters. 
colbar <- rainbow( max( memb$membership ) + 1 )
col <- colbar[ memb$membership + 1 ]
col[ is.na( col ) ] <- "grey"

lay <- layout.fruchterman.reingold( g, weights=E( g )$weight )
dev.new()
plot( g, layout=lay, vertex.size=3.5, vertex.label=V( g )$name, vertex.color=col,
		asp=FALSE, edge.color="gray", edge.width=1,
		main="hclust membership" )

## Output affils, but doesn't seem to capture weights.
#write.graph( g, "affils.net", format="ncol", weights=E(g)$weight )

## Copy names to ids for proper Pajek output.
V(g)$id <- V(g)$name
write.graph( g, "pajek.net", format="pajek" )
